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ABSTRACT 

The zero-temperature ordering kinetics of conserved XY models in spatial di- 
mensions d = 2 and 3 is studied using cell dynamical simulations. The growth of 
the characteristic length scale L(t) is fully consistent with recent theoretical predic- 
tions: L(t) ~ i 1//4 for d = 2 and L(t) ~ (tlnt) 1//4 for d — 3. A gaussian closure 
approximation describes the form of the structure factor rather well for d = 3. 
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There has been considerable recent interest in the kinetics of phase ordering in 
systems described by non-scalar order parameters The two fundamental ques- 
tions usually addressed are (a) the extent to which the familiar scaling phenomenol- 
ogy, developed in the context of scalar systems 0, is applicable to non-scalar sys- 
tems; and (b) if scaling holds, the nature of the growth law for the characteristic 
scale L{t) that describes the coarsening dynamics, t being the time. 

In recent work, Bray and Rutenberg (BR) have addressed the second of these 
questions using an 'energy scaling' argument ||. They obtain predictions for the 
asymptotic form of L(t) for all systems for which scaling holds and for which the 
dynamic is purely dissipative. 

The simplest non-scalar systems are XY models (or 0(2) systems). There have 
been a number of recent simulations of the ordering kinetics of XY systems, for both 
nonconserved 0, H and conserved [§, |7| dynamics, and some experimental studies 
of related liquid crystal systems ||. For nonconserved dynamics, the growth laws 
predicted by BR are L(t) ~ (t/lnt) 1 / 2 for d = 2, and L(t) ~ t 1/2 for d > 3 (see 
also related work by Yurke et al. ||), provided that scaling holds. For conserved 
dynamics, the corresponding predictions are |3j Lit) ~ t 1 / 4 for d = 2, and Lit) ~ 
(tint) 1 / 4 for d > 3. 

The present work has two goals. Firstly, we present numerical results from an 
extensive Cell Dynamical System (CDS) ]TD[ study of phase ordering dynamics for 
the conserved XY model in d = 2 and 3. In contrast to previous numerical works, 
we present data for the 'hardened structure factor', enabling a clear observation of 
the generalized Porod tail behavior |Tl|. Our results unambiguously demonstrate 
dynamical scaling for both d = 2 and d = 3, and are consistent with the theo- 
retical predictions for L(t). In particular, there is clear evidence for the predicted 
logarithmic correction to power-law growth for d = 3. 
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The second goal of this study is to examine the utility of gaussian closure schemes 
in computing the analytic form of the time-dependent structure factor. Our ana- 
lytic results using a gaussian closure scheme are in reasonable agreement with the 
numerical data for d = 3, but do not show good agreement with the data for d = 2. 
This suggests that, as is known to be true in the nonconserved case [H gaussian 
closure schemes for the conserved case are better for higher dimensionality 

The dynamical equation (in dimensionless form) for the ordering of the conserved 
XY model has the form 

dtf = -V 2 [V 2 - dV/d$] , (1) 

where is a 2-component vector order parameter and V(<j>) has a Mexican hat 
form, e.g. V(4>) = (1 — <f) 2 ) 2 , with a degenerate manifold of ground states, \<ft\ = 1. 
The CDS models used in our simulations were obtained by a conventional Euler 
discretization of (|l|) using an isotropic discrete Laplacian. The mesh sizes of our 
discretization were so large that it would be incorrect to claim that the solution 
of our numerical scheme accurately shadows that of the original partial differential 
equation. Thus, our models are justifiable only as CDS models belonging to the 
same dynamical universality class as the underlying partial differential equation 
|TIJ. Our d = 2 simulations were carried out on lattices of size 256 2 (with mesh 
sizes At = 0.15 and Ax = 1.7), and the spherically averaged structure factor was 
computed as an average over 80 runs with independent random initial conditions. 
Our d = 3 simulations were carried out on lattices of size 64 3 (with mesh sizes 
At = 0.1 and Ax = 1.7) and the structure factor was calculated as an average 
over 50 independent runs. The data were 'hardened', i.e. the structure factor was 
computed using fields renormalized to the length obtained from the fixed points of 
the CDS iteration scheme. This procedure gives better scaling at large momenta and 
elucidates the asymptotic tail behavior, which can be masked at late times by finite 
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defect sizes [jl3|] . If scaling is valid, the spherically averaged structure factor S(k,t) 
takes the scaling form S(k,t) = L d g(kL), where L = L(t) is the characteristic scale 
at time t, and g(x) is a universal scaling function. The length L(t) was defined in the 
usual way to be the reciprocal of the first moment, (k), of the structure factor, i.e. 
L(t) = (k) -1 , where (k) = mi/m , and the m n are the 'moments' of the structure 
factor: m n = J °° dk k n S(k, t). In this way, the length scale L{t) is extracted directly 
from the structure factor data. 

The scaled structure factor data are shown in Figures 1 (d = 3) and 2 (d = 2), 
in (a) linear-linear and (b) log-log form. The data collapse is good for both systems. 
The log-log plots reveal the large- and small-x behavior of the scaling function 
g(x). The large- x behavior follows the predicted power-law form [|T| g(x) ~ x -{d+2) 
('generalized Porod law'), which is a consequence of the vortex (d = 2) or vortex-line 
(d = 3) topological defects present in the system ||14|| . The continuous curves, which 
are obtained from the approximate analytical treatment discussed shortly, have this 
feature built in and thereby serve as useful guides to the eye in the large-x regime. 
There is a hump in the log-log plots of Figures 1(b) and 2(b) at x ~ 3. This is 
more clearly seen in the 'Porod plot' of g(x) x d+2 against x (which we will present 
elsewhere) and is reminiscent of the hump in the structure factor for the conserved 
scalar case ||15|| . For small x, the data are consistent with the expected x 4 behavior 
(broken lines), which can be derived |I| using an extension of the method used to 
obtain the x 4 behavior for scalar systems [16[ . 

The results for L[t) = (k) -1 are presented in Figure 3. In both cases (d = 2,3) 
we attempt two different fits, a simple power- law fit L(t) = At x (Figure 3(a)), and a 
power-law corrected by a logarithm, Lit) = A[t ln(t/r)] x (Figure 3(b)). The values 
of the best-fit exponent x in each case, and the timescale r for the logarithmic fits, 
are shown on the Figure. For d = 3, the logarithmic fit is extremely good, and much 
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better than a simple power-law. The best-fit exponent x = 0.250 is in perfect agree- 
ment with the BR prediction x = 1/4. When a simple power law is forced through 
the d = 3 data, the fit is much poorer than with the logarithmic correction, and the 
value of x (0.30) is unreasonably large. By contrast, for the d = 2 data a simple 
power-law gives a good fit (much better than for d — 3), with a best-fit exponent 
x = 0.247, again pleasingly close to the theoretical prediction x = 1/4. The loga- 
rithmic fit also works acceptably well for this case (as the presence of an additional 
fitting parameter, i.e. r, would lead us to expect), but the corresponding value of x 
(0.21) is unreasonably small. We conclude that the data for both dimensionalities 
are completely consistent with the BR prediction. 

Previous studies of the conserved XY model in d = 2 by Mondello and Goldenfeld 
H also found L(t) ~ t 25 , but there was some evidence for weak scaling violations: 
it was not possible to simultaneously collapse both the position and height of the 
peak in S(k,t). Of course, we should point out that the results presented in || 
had a substantially larger dynamic range (about 2.3 decades in t) than those in 
Figure 2 (one decade in t). However, the data in || exhibit weak scaling violations 
even over dynamic ranges comparable with ours. Apart from this, we note that our 
statistics (average over 80 runs) are somewhat better than those of Mondello and 
Goldenfeld (40 runs), and that their data were not hardened. However, is not clear 
to us that either of these facts should cause such an appreciable improvement around 
the peak position. In this context it is worth noting that recent studies p| of the 
nonconserved XY model in d — 2 found evidence for strong scaling violations when 
the characteristic spacing between vortices, d(t) = p -1 / 2 with p being the vortex 
density, was used as the scaling length L{t). It would be interesting to carry out 
a similar study for the conserved case (the vortex density was not measured in the 
present simulations). 
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Previous simulations for d = 3 conserved XY systems 0, on smaller systems 
(48 3 ) than those studied here, were originally interpreted in terms of power-law 
growth, L(t) ~ t x . For T = 0, the value x ~ 0.29 was obtained, close to the value 
x ~ 0.30 obtained here with a simple power-law fit, but the quality of the fit was not 
very good. Subsequently Siegert has shown that the logarithmic form proposed by 



BR (with exponent 1/4) gives a very good fit [17]. Some simulations at T > were 
fitted to a power-law with x ~ 0.26, but again the logarithmic form (with exponent 
1/4) gives a much better fit (but with a different r than for T = 0) |T7| . 

We turn now to our approximate analytical treatment, based on a 'gaussian 
closure approximation' applied to the equation of motion. This approach requires a 
straightforward modification of equivalent treatments of nonconserved n- vector fields 



fli8| , which in turn generalize an earlier treatment of nonconserved scalar fields |L9 
Here we just sketch the derivation and obtain results for the conserved XY model 
(n = 2). A fuller treatment, with applications to other values of n, will be given in a 



longer publication ||20|| . The starting point is the extension ([!]) of the Cahn-Hilliard 



equation (T[| to vector fields. The first step is to take the scalar product of (0), 
evaluated at space-time point '1', with 0(2), the field at space-time point '2', and 
average over an ensemble of initial conditions. This gives 

d tl C(12) = -V 2 [V 2 C(12) - (0(2) • dV/d${l))\ , (2) 

where C(12) = (0(1) • 0(2)) is the pair correlation function. In deriving (0), we 
exploited the translational invariance of the ensemble of initial conditions. Following 
earlier studies of nonconserved systems ||18|| , we impose an approximate closure of 



the exact equation (||) through a gaussian assumption for an auxiliary vector field 
m(x, t), related to the physical field via through the equation 

V 2 J =dV/d$ (3) 



for the function c/>(m). Eq. @ is to be solved with boundary conditions 0(0) = 0, 
<f)(m) — * rh for \rh\ — > oo, where m = m/|m| is a unit vector. The Laplacian operator 
in Eq. (|3]) is defined by = Yh=i d 2 /dm 2 , where the m^'s are the Cartesian 
components of the vector rh. The physical meaning of (Q) is that the function c/>(m) 
gives the structure of an equilibrium defect, with \fh\ regarded as the distance from 



the defect core |18 



The key approximation here is to treat rh as a gaussian random field in (Q). 
For nonconserved fields, this approximation is qualitatively accurate, and becomes 
quantitatively accurate with increasing d |J], |T2f . For conserved fields its status 
is less clear, but it provides a useful starting point for the discussion of scaling 
functions for conserved vector fields. (We should warn the reader that the gaussian 
approximation in its present form is not even a reasonable starting point for the 
conserved scalar case, which has an altogether different growth law, i.e. L{t) ~ t 1 / 3 , 
from the conserved vector case.) Using (Q) to eliminate the explicit dependence 
on the potential in (|), and exploiting the gaussian property, we obtain a closed 
equation for C(12). For equal-time correlations, this reads 



ldC ^2 
= -V 2 

2 dt 



1 d C 

{m z ) cry 



(4) 

In (01) (m 2 ) is the variance of one component m of rh, 7 = (m(l)m(2)) / (m 2 ) is the 



normalized correlator of m, and the function C(j) for general n is given by fTI 
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where -B(x, ?/) is the beta function and 2 Fi (a, 6; c; z) is the Hypergeometric function. 

The next step is to seek a scaling solution C = f(r/L), in which all the time 
dependence is contained in L = L(t). Requiring that all terms in ([|) scale in the 
same way forces (m 2 ) = aL 2 and L = (St) 1 / 4 (where the factor 8 is put in for 
convenience). Recasting (^) as an equation for 7(x), where x = r/L is the scaling 



variable, gives 



|v z - 2 (/ 7 x 7 ') = vh + « 7 + f -f {if , (6) 



which is a convenient form for numerical solution. Here the function 7(7) is just the 
right-hand side of (||), subscripts 7 indicate derivatives with respect to 7, primes 
denote derivatives with respect to x, and V 2 X is the Laplacian with respect to x, i.e. 
V 2 X = d 2 /dx 2 + \{d - l)/x]d/dx. 

In an earlier treatment of Eqs. (Q) and (Q), an approximate solution was obtained 
by expanding •ydC/d'y in powers of C and truncating at 0(C 3 ) |^T|. The Porod tail 
in the structure factor is lost in this approximation. Despite this, the resulting real- 
space correlation function was in good agreement with the data of Siegert and Rao 
0, although this agreement may be misleading as the data in [0] were not hardened. 
The present treatment is superior, as no approximations are made beyond the initial 
gaussian closure. In particular, the Porod tail is present in the solution, as is evident 
from Figures 1(b) and 2(b). 

The left-hand side of (BJ) can be written in integral form using the identity 

1 



V-F(x) 



d-2 



J o X dyyF(y)U^j d * - lj + j™ dyyF(y) 



(7) 



for d > 2. (For d = 2 the appropriate limit has to taken p(j). In the context of Eq. 
(D,F(x)=/ 7 xf(x)=x/ / (x) in®. 

To solve Eq. (§) numerically requires specifying boundary conditions at x = 0. 
In this paper we specialise to the case n = 2. Solutions for other values of n will be 
presented elsewhere [p0| . For n = 2 one can show that 7(2;) has a small-x expansion 
of the form 

tf*> = i-5«V-£; + - < 8 > 

Therefore the required boundary conditions are 7(0) = 1, 7'(0) = 0. To perform 
the numerical integration, however, one needs the values of the parameter a in (||) 



and of the infinite integral J °° dyyF(y) in (0). Using F(x) = xf'(x), this integral 
can be written as —2 f£°dyyf(y), which can be expressed in terms of a and the 
parameters a and b in the small- x expansion @. Inserting (§) into ([]), using (^) 
to simplify the left-hand side, and matching terms of 0(1) and 0(1/ Inx) gives the 
relations 

a = da 2 (9) 

POO 

2 / dyyf(y) = (d - 2)(a 2 + 2bd) . (10) 
Jo 

(Again, the case d = 2 has to be treated separately |2"0]]). Numerical integration 
of the differential equation @ is now straightforward, as the finite integral in (0) 
can be accumulated as the integration proceeds. The two unknown parameters a 
and b are fixed by the requirement that the large- a; behavior be physical, i.e. f(x) 
must vanish as x — > oo. For d = 3, these parameters take the approximate values 
a = 1.177567, b = —0.155217. (For d = 2, there are still two parameters, but they 
enter in a slightly different way |20[ ) . 

The numerical solutions for the structure factor (the Fourier transform of the 
pair correlation function) are included in Figures 1 and 2 as solid lines. Notice that 
there are no adjustable parameters in the fit. The fit is surprisingly good for d = 3, 
despite the fact that the growth law L ~ t 1 / 4 obtained within the gaussian closure 
scheme does not have the logarithmic correction predicted for the exact solution. 
For d = 2, the fit is noticeably worse. One feature of the approximate theory which 
is qualitatively incorrect is the small- k behavior. The log-log plots of Figs. 2(b) and 
3(b) show the expected k A behavior at small k. This feature is missing from the 
gaussian closure scheme, which gives a k 2 behavior at small k. Another point of 
discrepancy between the gaussian closure theory and the numerical results is that 
the theory overestimates the amplitude of the Porod tail though it replicates nicely 
the qualitative features (i.e. valleys and humps), especially in d = 3. This feature is 
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better seen in a Porod plot of theory and data, and we will present such a plot in 
an extended publication. 

In summary, our simulation data for the conserved XY model is fully consistent, 
for both d = 2 and d = 3, with theoretical predictions derived from a scaling 
assumption ||. The scaling collapse of the data, using a lengthscale (/c) -1 extracted 
from the structure factor data itself, is good. In future work it would be desirable to 
test scaling using a lengthscale derived independently from, for example, the density 
of vortices (or vortex lines) ||. An approximate calculation of the structure factor 
scaling function, using a gaussian closure scheme, describes the data quite well for 
d = 3, except at small k. We expect this procedure to work even better for n > 2, 
where the growth law L ~ i 1 / 4 has no logarithmic corrections || . 

SP thanks the members of the Physics Department of Manchester University, 
where most of this work was done, and he and AB thank the staff of the Isaac 
Newton Institute, Cambridge, where the work was completed, for their hospitality. 
FR thanks CONACYT (Mexico) for support. 



9 



References 

[1] For a recent review see A. J. Bray, Adv. Phys. 43, 357 (1994). 

[2] K. Binder and D. Stauffer, Phys. Rev. Lett. 33, 1006 (1974); J. Marro, J. L. 
Lebowitz and M. H. Kalos, Phys. Rev. Lett. 43, 282 (1979); H. Furukawa, Prog. 
Theor. Phys. 59, 1072 (1978); Phys. Rev. Lett. 43, 136 (1979). 

[3] A. J. Bray and A. D. Rutenberg, Phys. Rev. E 49, R 27 (1994); A. D. Rutenberg 
and A. J. Bray, Phys. Rev. E, in press. 

[4] M. Mondello and N. Goldenfeld, Phys. Rev. A 42, 5865 (1990); Phys. Rev. A 
45, 657 (1992); H. Toyoki, J. Phys. Soc. Jpn. 60, 1433 (1991). 

[5] R. E. Blundell and A. J. Bray, Phys. Rev. E 49, 4925 (1994). 

[6] M. Mondello and N. Goldenfeld, Phys. Rev. E 47, 2384 (1993). 

[7] M. Siegert and M. Rao, Phys. Rev. Lett. 70, 1956 (1993). 

[8] A. N. Pargellis, S. Green and B. Yurke, Phys. Rev. E 49, 4250 (1994). 

[9] B. Yurke, A. N. Pargellis, T. Kovacs and D. A. Huse, Phys. Rev. E 47, 1525 
(1993). 

[10] Y. Oono and S. Puri, Phys. Rev. Lett. 58, 836 (1987); Phys. Rev. A 38, 434 
(1988); S. Puri and Y. Oono, Phys. Rev. A 38, 1542 (1988). 

[11] A. J. Bray and S. Puri, Phys. Rev. Lett. 67, 2670 (1991); H. Toyoki, Phys. 
Rev. B 45, 1965 (1992). 

[12] A. J. Bray and K. Humayun, Phys. Rev. E 48, R1609 (1993). 

[13] Y. Oono and S. Puri, Mod. Phys. Lett. B 2, 861 (1988). 

10 



[14] A. J. Bray and K. Humayun, Phys. Rev. E 47, R9 (1993). 

[15] A. Shinozaki and Y. Oono, Phys. Rev. Lett. 66, 173 (1991). 

[16] C. Yeung, Phys. Rev. Lett. 61, 1135 (1988); H. Furukawa, Phys. Rev. B 40, 
2341 (1989); Prog. Theor. Phys. Suppl. 99, 358 (1989); H. Tomita, Prog. Theor. 
Phys. 85, 47 (1991); P. Fratzl, J. L. Lebowitz, O. Penrose and J. Amar, Phys. 
Rev. B 44, 4794 (1991). 

[17] M. Siegert, to appear in Scale Invariance, Interfaces, and N on- Equilibrium 
Dynamics, proceedings of the NATO Advanced Study Institute, edited by A. 
J. McKane, M. Droz, J. Vannimenus and D. Wolf (Plenum, New York, 1994). 

[18] Fong Liu and G. F. Mazenko, Phys. Rev. B 45, 6989 (1992); A. J. Bray and K. 
Humayun, J. Phys. A 25, 2191 (1992). 

[19] G. F. Mazenko, Phys. Rev. Lett. 63, 1605 (1989); Phys. Rev. B 42, 4487 (1990). 

[20] S. Puri, A. J. Bray and F. Rojas, unpublished. 

[21] F. Rojas and A. J. Bray, Phys. Rev. E 51, 188 (1995). 

[22] G. F. Mazenko, University of Chicago preprint. 



11 



Figure Captions 

1. Scaled structure factor for d = 3: (a) linear-linear plot (b) log-log plot. The 
data are for dimensionless times 1000, 2000, 3000 and 5000, as indicated. 
Details of the simulation are described in the text. The continuous curves are 
the numerical Fourier transforms of the scaling function for the pair correlation 
function, calculated using the gaussian closure scheme described in the text. 
The dashed line in (b) has slope 4. 

2. Same as Figure 2, but for d = 2. Here the data are for dimensionless times 
1050, 2100, 6300 and 10500. 

3. Growth of the characteristic length scale L{t) for the conserved XY model, 
fitting to (a) L{t) = At x and (b) L(t) = A(t \n[t/r]) x . The best-fit values of x 
and r are specified on the figures. 



12 



